package scu.maqiang.mesh;

import static java.lang.Math.abs;

import java.io.*;
import java.lang.reflect.Array;
import java.util.*;
import java.util.function.Function;
import java.util.stream.IntStream;
import java.util.stream.Stream;

import scu.maqiang.numeric.*;

/**
 * Mesh File
 * 
 * Dimension
 * 2
 * 
 * NumberOfNodesInOneElement
 * 3
 * 
 * NumberOfNodesOnBoundary
 * 2
 * 
 * Nodes
 * 15
 * x1 y1 lbl1
 * x2 y2 lbl2
 * ...
 * x15 y15 lbl15
 * 
 * Elements
 * 20
 * e11 e12 e13 lbl1
 * e21 e22 e23 lbl2
 * ...
 * e201 e202 e203 lbl20
 * 
 * Boundarys
 * 8
 * b11 b12 lbl1
 * b21 b22 lbl2
 * ...
 * b81 b82 lbl8
 * 
 * End
 */

/**
 * 有限元抽象网格类，定义有限元中所有网格都应该具有的一些属性与方法
 * 
 * @author 马强
 *
 */
public abstract class Mesh {

	public Mesh() {

	}
	/**
	 * 得到网格结点个数
	 * 
	 * @return 网格结点个数
	 */
	public int getNv() {
		return nv;
	}

	/**
	 * 得到网格单元个数
	 * 
	 * @return 网格单元个数
	 */
	public int getNt() {
		return nt;
	}

	public abstract int getTecplotNt();

	/**
	 * 得到网格维度
	 * 
	 * @return 网格维数
	 */
	public int getnDim() {
		return nDim;
	}

	/**
	 * 初始化网格
	 * 
	 * @param nv 网格结点个数
	 * @param nt 网格区域单元个数
	 * @param nb 网格边界单元个数
	 */
	public void initMesh0(int nv, int nt, int nb) {
		this.nv = nv;
		this.nt = nt;
		this.nb = nb;
		this.nodes = new double[nv][];
		this.nodeLabel = new int[nv];
		this.elements = new int[nt][];
		this.elementLabel = new int[nt];
		this.boundaries = new int[nb][];
		this.boundaryLabel = new int[nb];
		this.comparators = new Comparator[3];
		comparators[0] = this.comparator1;
		comparators[1] = this.comparator2;
		comparators[2] = this.comparator3;
	}

	/**
	 * 得到网格上第个单元的结点个数
	 * 
	 * @return 每个单元的结点个数
	 */
	public int getnPerElement() {
		return nPerElement;
	}

	/**
	 * 得到网格每个边界单元的结点个数
	 * 
	 * @return 边界单元结点个数
	 */
	public int getnPerBoundary() {
		return nPerBoundary;
	}

	/**
	 * 得到网格边界单元个数
	 * 
	 * @return 网格边界单元个数
	 */
	public int getNb() {
		return nb;
	}

	/**
	 * 得到网格对应的tecplot网格类型名
	 * 
	 * @return ecplot网格类型名
	 */
	public String getTecplotType() {
		return tecplotType;
	}

	/**
	 * 得到网格结点流
	 * 
	 * @return
	 */
	public Stream<double[]> getNodeStream() {
		return Stream.of(nodes);
	}

	/**
	 * 得到网格区域单元流
	 * 
	 * @return
	 */
	public Stream<int[]> getElementStream() {
		return Stream.of(elements);
	}

	/**
	 * 得到网格边界单元流
	 * 
	 * @return
	 */
	public Stream<int[]> getBoundaryStream() {
		return Stream.of(boundaries);
	}

	/**
	 * 计算网格中单元结点编号最大差值，用于计算形成的稀疏矩阵半带宽
	 * 
	 * @return 单元结点编号最大差值
	 */
	public int elementDiff() {
		int result = 0;
		for (int k = 0; k < nt; k++) {
			int[] idx = elements[k];
			int n = idx.length;
			for (int i = 0; i < n; i++) {
				for (int j = i + 1; j < n; j++) {
					int diff = Math.abs(idx[i] - idx[j]);
					if (result < diff) {
						result = diff;
					}
				}
			}
		}
		return result;
	}

	public void constructMesh(double[][] coords, int[][] eles, boolean zeroBased) {
		assert (coords[0].length >= nDim);
		int base = zeroBased? 0: -1;
		nv = coords.length;
		nodes = new double[nv][nDim];
		nodeLabel = new int[nv];
		for(int i = 0; i < nv; i++) {
			for(int j = 0; j < nDim; j++) {
				nodes[i][j] = coords[i][j];
			}
		}

		nt = eles.length;
		elements = new int[nt][nPerElement];
		elementLabel = new int[nt];
		for(int i = 0; i < nt; i++) {
			for(int j = 0; j < nPerElement; j++) {
				elements[i][j] = eles[i][j] + base;
			}
		}
		extractBoundarysFromElements();
	}
//	public void readMesh(MeshFormat meshType, String fileName) {
//		switch (meshType) {
//			case msh -> readMesh_Msh(fileName);
//			case mesh -> readMesh_Mesh(fileName);
//			default -> throw new IllegalArgumentException("Unsupported MeshFormat");
//		}
//
//	}

//	private void readMesh_Msh(String fileName) {
//		Scanner sc = null;
//		String line = "";
//		int nv, nt, nb;
//		try {
//			sc = new Scanner(new File(fileName));
//		} catch (FileNotFoundException e) {
//			e.printStackTrace();
//		}
//		nv = sc.nextInt();
//		nt = sc.nextInt();
//		nb = sc.nextInt();
//		initMesh0(nv, nt, nb);
//		for (int i = 0; i < nv; i++) {
//			double x = sc.nextDouble();
//			double y = sc.nextDouble();
//			nodes[i] = new double[] {x, y};
//			nodeLabel[i] = sc.nextInt();
//		}
//		for (int i = 0; i < nt; i++) {
//			elements[i] = new int[nPerElement];
//			for (int j = 0; j < nPerElement; j++) {
//				elements[i][j] = sc.nextInt() - 1;
//			}
//			elementLabel[i] = sc.nextInt();
//		}
//		for(int i = 0; i < nb; i++) {
//			boundaries[i] = new int[nPerBoundary];
//			for(int j = 0; j < nPerBoundary; j++) {
//				boundaries[i][j] = sc.nextInt() - 1;
//			}
//			boundaryLabel[i] = sc.nextInt();
//		}
//		sc.close();
//	}

//	private void readMesh_Mesh(String fileName) {
//		try (BufferedReader br = new BufferedReader(new FileReader(fileName))) {
//			String ll = br.readLine();
//			while (!"Vertices".equals(ll.trim())) {
//				ll = br.readLine();
//			}
//			nv = Integer.parseInt(br.readLine());
//			nodes = new double[nv][];
//			nodeLabel = new int[nv];
//			double[] coord = new double[3];
//			for (int i = 0; i < nv; i++) {
//				String[] data = br.readLine().split("\\s+");
//				nodes[i] = new double[nDim];
//				for(int j = 0; j < nDim; j++) {
//					nodes[i][j] = Double.parseDouble(data[j]);
//				}
//				nodeLabel[i] = Integer.parseInt(data[3]);;
//			}
//			while (!meditMeshType[1].equals(ll.trim())) {
//				ll = br.readLine();
//			}
//			nt = Integer.parseInt(br.readLine());
//			elements = new int[nt][];
//			elementLabel = new int[nt];
//			for (int i = 0; i < nt; i++) {
//				String[] data = br.readLine().split("\\s+");
//				elements[i] = new int[nPerElement];
//				for(int j = 0; j < nPerElement; j++) {
//					elements[i][j] = Integer.parseInt(data[j]) - 1;
//				}
//				elementLabel[i] = Integer.parseInt(data[nPerElement]);;
//			}
//			while (!meditMeshType[2].equals(ll)) {
//				ll = br.readLine();
//			}
//			nb = Integer.parseInt(br.readLine());
//			nb = nb;
//			boundaries = new int[nb][];
//			boundaryLabel = new int[nb];
//			for (int i = 0; i < nb; i++) {
//				String[] data = br.readLine().split("\\s+");
//				boundaries[i] = new int[nPerBoundary];
//				for(int j = 0; j < nPerBoundary; j++) {
//					boundaries[i][j] = Integer.parseInt(data[j]) - 1;
//				}
//				boundaryLabel[i] = Integer.parseInt(data[nPerBoundary]);
//			}
//		} catch (IOException e) {
//			e.printStackTrace();
//		}
//	}

	public void readMesh(MeshFormat meshType, String fileName) {
		switch (meshType) {
			case msh -> readMesh_Msh(fileName);
			case mesh -> readMesh_Mesh(fileName);
			default -> throw new IllegalArgumentException("Unsupported MeshFormat");
		}
	}

	private void readMesh_Msh(String fileName) {
		Scanner sc = null;
		String line = "";
		int nv, nt, nb;
		try {
			sc = new Scanner(new File(fileName));
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
		nv = sc.nextInt();
		nt = sc.nextInt();
		nb = sc.nextInt();
		initMesh0(nv, nt, nb);
		for (int i = 0; i < nv; i++) {
			nodes[i] = new double[nDim];
			for(int j = 0; j < nDim; ++j) {
				nodes[i][j] = sc.nextDouble();
			}
			nodeLabel[i] = sc.nextInt();
		}
		for (int i = 0; i < nt; i++) {
			elements[i] = new int[nPerElement];
			for(int j = 0; j < nPerElement; ++j) {
				elements[i][j] = sc.nextInt() - 1;
			}
			elementLabel[i] = sc.nextInt();
		}
		for(int i = 0; i < nb; ++i) {
			boundaries[i] = new int[nPerBoundary];
			for(int j = 0; j < nPerBoundary; ++j) {
				boundaries[i][j] = sc.nextInt() - 1;
			}
			boundaryLabel[i] = sc.nextInt();
		}
		sc.close();
	}

	private void readMesh_Mesh(String fileName) {
		try (BufferedReader br = new BufferedReader(new FileReader(fileName))) {
			String ll = br.readLine();
			while (!meditMeshType[0].equals(ll.trim())) {
				ll = br.readLine();
			}
			nv = Integer.parseInt(br.readLine());
			nodes = new double[nv][];
			nodeLabel = new int[nv];
			double[] coord = new double[3];
			for (int i = 0; i < nv; i++) {
				String[] data = br.readLine().split("\\s+");
				coord[0] = Double.parseDouble(data[0]);
				coord[1] = Double.parseDouble(data[1]);
				coord[2] = Double.parseDouble(data[2]);
				int num = Integer.parseInt(data[3]);
				nodes[i] = new double[nDim];
				for(int j = 0; j < nDim; j++) {
					nodes[i][j] = coord[j];
				}
				nodeLabel[i] = num;
			}
			while (!meditMeshType[1].equals(ll.trim())) {
				ll = br.readLine();
			}
			nt = Integer.parseInt(br.readLine());
			elements = new int[nt][];
			elementLabel = new int[nt];
			int[] ele = new int[nPerElement];
			for (int i = 0; i < nt; i++) {
				String[] data = br.readLine().split("\\s+");
				for(int j = 0; j < nPerElement; ++j) {
					ele[j] = Integer.parseInt(data[j]) - 1;
				}
				elements[i] = ele.clone();
				elementLabel[i] = Integer.parseInt(data[nPerElement]);
			}
			while (!meditMeshType[2].equals(ll)) {
				ll = br.readLine();
			}
			nb = Integer.parseInt(br.readLine());
			boundaries = new int[nb][];
			boundaryLabel = new int[nb];
			ele = new int[nPerBoundary];
			for (int i = 0; i < nb; i++) {
				String[] data = br.readLine().split("\\s+");
				for(int j = 0; j < nPerBoundary; ++j) {
					ele[j] = Integer.parseInt(data[j]) - 1;
				}
				boundaries[i] = ele.clone();
				boundaryLabel[i] = Integer.parseInt(data[nPerBoundary]);
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	public void saveMesh(MeshFormat meshType, String fileName) {
		switch (meshType) {
			case msh -> saveMesh_Msh(fileName);
			case mesh -> saveMesh_Mesh(fileName);
			default -> throw new IllegalArgumentException("Unsupported MeshFormat");
		}
	}
	/**
	 * 将网格文件保存到文件中，网格文件后缀名为.msh
	 * 
	 * @param fileName 保存的网格文件名
	 */
	private void saveMesh_Msh(String fileName) {
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
			bw.write(nv + "\t" + nt + "\t" + nb + "\n");
			for (int i = 0; i < nv; i++) {
				double[] node = nodes[i];
				for (int j = 0, n = node.length; j < n; j++) {
					bw.write(node[j] + "\t ");
				}
				bw.write(nodeLabel[i] + "\n");
			}

			for (int i = 0; i < nt; i++) {
				int[] ele = elements[i];
				for (int j = 0, n = ele.length; j < n; j++) {
					bw.write((ele[j] + 1) + "\t ");
				}
				bw.write(elementLabel[i] + "\n");
			}

			for (int i = 0; i < nb; i++) {
				int[] ele = boundaries[i];
				for (int j = 0, n = ele.length; j < n; j++) {
					bw.write((ele[j] + 1) + "\t ");
				}
				bw.write(boundaryLabel[i] + "\n");
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	private void saveMesh_Mesh(String fileName) {
		try(BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
			bw.write("MeshVersionFormatted 1\n\n");
			bw.write("Dimension " + nDim + "\n");
			bw.write(meditMeshType[0] + "\n");
			bw.write(nv + "\n");
			for(int i = 0; i < nv; ++i) {
				double[] node = nodes[i];
				for(int j = 0, nj = node.length; j < nj; ++j) {
					bw.write(node[j] + "  ");
				}
				for(int j = node.length; j < nDim; j++) {
					bw.write("0.0  ");
				}
				bw.write(nodeLabel[i] + "\n");
			}
			bw.write("\n");
			bw.write(meditMeshType[1] + "\n");
			bw.write(nt + "\n");
			for(int i = 0; i < nt; i++) {
				int[] ele = elements[i];
				for(int j = 0, nj = ele.length; j < nj; j++) {
					bw.write((ele[j] + 1) +"  ");
				}
				bw.write(elementLabel[i] + "\n");
			}
			bw.write("\n");
			bw.write(meditMeshType[2] + "\n");
			bw.write(nb + "\n");
			for(int i = 0; i < nb; i++) {
				int[] ele = boundaries[i];
				for(int j = 0, nj = ele.length; j < nj; j++) {
					bw.write((ele[j] + 1) + "  ");
				}
				bw.write(boundaryLabel[i] + "\n");
			}
			bw.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	/**
	 * 设置节点材料标识
	 * @param i 节点编号
	 * @param label 新的节点标识
	 */
	public void setNodeLabel(int i, int label) {
		nodeLabel[i] = label;
	}

	/**
	 * 设置节点材料标识
	 */
	public void setNodeLabel(ScalarFunc regionFunc, double[] param, int label) {
		for (int i = 0; i < nt; i++) {
			int[] ele = elements[i];
			if (regionFunc.action(interiorPointInElement(ele), elementLabel[i], param) > 0) {
				for(int k = 0, n = ele.length; k < n; k++) {
					nodeLabel[ele[k]] = label;
				}
			}
		}
	}

	/**
	 * 设置结点坐标
	 * @param direct 设置方向维度
	 * @param u 新的坐标值
	 */
	public void setNodeCoordinates(Direct direct, double[] u) {
		ParamCheck.checkEqual(u.length, nv);
		int directVal = direct.getValue();
		ParamCheck.checkBounds(directVal, 0, nDim);
		for(int i = 0; i < nv; i++) {
			nodes[i][directVal] = u[i];
		}
	}

	public int getNodeLabel(int i) {
		return nodeLabel[i];
	}

	public int getDomainLabel(double[] xyz) {
		PointInElement pie = ConstantMaps.elementPointJudgeMap.get(this.getClass().getName());
		double[][] coord = new double[nPerElement][nDim];
		double[] N = new double[nPerElement];
		for(int i = 0; i < nt; i++) {
			getCoordPerElement(elements[i], coord);
			if(pie.action(coord, xyz, N)) {
				return elementLabel[i];
			}
		}
		throw new IllegalArgumentException("Point is not in domain!");
	}

	
	public double[] functionInMesh(ScalarFunc func, double[] param) {
		double[] u = new double[nv];
		double[] coordNode = new double[nDim];
		for (int i = 0, n = u.length; i < n; i++) {
			getCoordPerNode(i, coordNode);
			u[i] = func.action(coordNode, 0, param);
		}
		return u;
	}

	public Complex[] functionInMesh(ScalarCFunc func, Complex[] param) {
		Complex[] u = new Complex[nv];
		double[] coordNode = new double[nDim];
		for (int i = 0, n = u.length; i < n; i++) {
			getCoordPerNode(i, coordNode);
			u[i] = func.action(coordNode, 0, param);
		}
		return u;
	}

	public void setFuncValueInNode(ScalarFunc func, double[] param, double value, double[] u) {
		double[] coordNode = new double[nDim];
		for (int i = 0, n = u.length; i < n; i++) {
			getCoordPerNode(i, coordNode);
			if(func.action(coordNode, nodeLabel[i], param) > 0.0) {
				u[i] = value;
			}
		}
	}
	public void setFuncValueInElement(ScalarFunc func, double[] param, double value, double[] u) {
		double[][] coordElement=  new double[nPerElement][nDim];
		double[] point;
		int[] ele;
		for (int i = 0; i < nt; i++) {
			ele = elements[i];
			getCoordPerElement(ele, coordElement);
			point = getinteriorPointInElement(coordElement);
			if(func.action(point, elementLabel[i], param) > 0.0) {
				for(int j = 0; j < nPerElement; j++) {
					u[ele[j]] = value;
				}
			}			
		}
	}

	public void setValueOnBoundary(double[] u, double val, int... bdLabel) {
		HashSet<Integer> bdSet = extractBoundaryNodes(bdLabel);
		for(int idx: bdSet) {
			u[idx] = val;
		}
	}

	public void setValueOnBoundary(ScalarFunc func, double[] param, double[] u, int... bdLabel) {
		HashSet<Integer> bdSet = extractBoundaryNodes(bdLabel);
		double[] coord = new double[nDim];
		for(int idx: bdSet) {
			getCoordPerNode(idx, coord);
			u[idx] = func.action(coord, 0, param);
		}
	}

	public void setFuncValueOnBoundary(ScalarFunc func, double[] param, double value, double[] u) {
		double[][] coordBoundary = new double[nPerBoundary][nDim];
		double[] point;
		int[] ele;
		for (int i = 0; i < nt; i++) {
			ele = boundaries[i];
			getCoordPerBoundary(ele, coordBoundary);
			point = getinteriorPointOnBoundary(coordBoundary);
			if (func.action(point, boundaryLabel[i], param) > 0) {
				for(int j = 0; j < nPerBoundary; j++) {
					u[ele[j]] = value;
				}	
			}		
		}
	}
	
	
	public boolean isPointInDomain(double[] xyz, double[] N) {
		PointInElement pie = getPointJudgeFunc();
		double[][] coord = new double[nPerElement][nDim];
		for(int i = 0; i < nt; i++) {
			getCoordPerElement(elements[i], coord);
			if(pie.action(coord, xyz, N)) {
				return true;
			}
		}
		return false;
	}
	
	public PointInElement getPointJudgeFunc() {
		return ConstantMaps.elementPointJudgeMap.get(this.getClass().getName());
	}
	
	public void setDomainLabel(int i, int label) {
		elementLabel[i] = label;
	}

	public void addDomainAsBoundaries(int domainLabel, int newBdLabel) {
		var edges = extractEdgesFromElements(domainLabel);
		addBoundarys(edges, newBdLabel);
	}

	public void addDomainAsBoundaries(ScalarFunc domainFunc, double[] param, int newBdLabel) {
		var edges = extractEdgesFromElements(domainFunc, param);
		addBoundarys(edges, newBdLabel);
	}

	/**
	 * 设置区域单元编号，用于定义各区域的材料编号
	 * 
	 * @param func  区域函数，对于要设置的单元，返回正数，否则返回负数
	 * @param param 区域函数参数
	 * @param label 要设置的单元区域编号
	 */
	public void setDomainLabel(ScalarFunc func, double[] param, int label) {
		for (int i = 0; i < nt; i++) {
			int[] ele = elements[i];
			if (func.action(interiorPointInElement(ele), elementLabel[i], param) > 0) {
				elementLabel[i] = label;
			}
		}
	}
	
	public void setBoundaryLabel(int i, int label) {
		boundaryLabel[i] = label;
	}

	/**
	 * 设置边界单元编号，用于区分边界以便设置边界条件
	 * 
	 * @param func  区域函数, 对于要设置的边界，返回正数，否则返回负数
	 * @param param 区域函数参数
	 * @param label 要设置的边界区域编号
	 */
	public void setBoundaryLabel(ScalarFunc func, double[] param, int label) {
		for (int i = 0; i < nb; i++) {
			int[] ele = boundaries[i];
			if (func.action(interiorPointOnBoundary(ele), boundaryLabel[i], param) > 0) {
				boundaryLabel[i] = label;
				//System.out.println("OK");
			}
		}
	}

	/**
	 * 得到单元中某一内点， 一般得到单元几何中心
	 * 
	 * @param ele 单元编号
	 * @return 第i个单元内一点坐标
	 */
	public double[] interiorPointInElement(int[] ele) {
		double[][] coord = new double[nPerElement][nDim];
		getCoordPerElement(ele, coord);
		return getinteriorPointInElement(coord);
	}

	public double[][] getCoords() {
		return MVO.transpose(nodes);
	}

	/**
	 * 得到边界中某一内点， 一般得到边界几何中心
	 * 
	 * @param ele 边界编号
	 * @return 第i个边界内一点坐标
	 */
	public double[] interiorPointOnBoundary(int[] ele) {
		double[][] coord = new double[nPerBoundary][nDim];
		getCoordPerBoundary(ele, coord);
		return getinteriorPointOnBoundary(coord);
	}

	/**
	 * 在得到单元各结点坐标基础之上，计算单元内某一点
	 * 
	 * @param coord 单元结点坐标矩阵，行数为单元结点个数nPerElement, 列数为网格维数nDim
	 * @return 单元内部点坐标，维数为nDim
	 */
	public abstract double[] getinteriorPointInElement(double[][] coord);

	/**
	 * 在得到边界各结点坐标基础之上，计算边界上某一点
	 * 
	 * @param coord 边界结点坐标矩阵，行数为单元结点个数nPerBoundary, 列数为网格维数nDim
	 * @return 边界上点坐标，维数为nDim
	 */
	public abstract double[] getinteriorPointOnBoundary(double[][] coord);

	/**
	 * 得到单元对应的边界信息
	 * 
	 * @param i 单元编号
	 * @return 锟竭斤拷锟斤拷锟斤拷锟斤拷锟斤拷
	 */
	public abstract int[][] getBoundarysFromElement(int i);

	/**
	 * 得到网格中每个单元相邻的单元列表, 使用一个HashMap记录网格各边界对应的单元编号,
	 *
	 * @return
	 */
	public int[][] getNeighbours() {
		int[][] neighbours = new int[nt][nBoundaryPerElement];
		for(int i = 0; i < nt; i++) {
			Arrays.fill(neighbours[i], -1);
		}
		HashMap<UnorderIntArray, Integer> elementRecord= new HashMap<>();
		for(int i = 0; i < nt; i++) {
			int[][] bds = getBoundarysFromElement(i);
			for(int j = 0; j < nBoundaryPerElement; j++) {
				UnorderIntArray ay = new UnorderIntArray(bds[j]);
				if (elementRecord.containsKey(ay)) {
					int eleInfo = elementRecord.get(ay);
					int ele = eleInfo / nBoundaryPerElement;
					int idx = eleInfo - nBoundaryPerElement * ele;
					neighbours[ele][idx] = i;
					neighbours[i][j] = ele;
					elementRecord.remove(ay);
				} else {
					elementRecord.put(ay, i * nBoundaryPerElement + j);
				}
			}
		}
		return neighbours;
	}

	/**
	 * 打印结点信息
	 */
	public void displayNodes() {
		for (int i = 0; i < nv; i++) {
			double[] node = nodes[i];
			for (int j = 0, n = node.length; j < n; j++) {
				System.out.print(node[j] + "\t ");
			}
			System.out.print(nodeLabel[i] + "\n");
		}
	}

	/**
	 * 打印区域单元信息
	 */
	public void displayElements() {
		for (int i = 0; i < nt; i++) {
			int[] ele = elements[i];
			for (int j = 0, n = ele.length; j < n; j++) {
				System.out.print(ele[j] + "\t ");
			}
			System.out.print(elementLabel[i] + "\n");
		}
	}

	/**
	 * 打印边界单元信息
	 */
	public void displayBoundarys() {
		for (int i = 0; i < nb; i++) {
			int[] ele = boundaries[i];
			for (int j = 0, n = ele.length; j < n; j++) {
				System.out.print(ele[j] + "\t");
			}
			System.out.print(boundaryLabel[i] + "\n");
		}
	}

	/**
	 * 得到第i个结点的坐标
	 * 
	 * @param i     结点编号
	 * @param coord 存储结点坐标
	 */
	public void getCoordPerNode(int i, double[] coord) {
		double[] temp = nodes[i];
		ParamCheck.checkEqual(nDim, coord.length);
		for (int j = 0; j < nDim; j++) {
			coord[j] = temp[j];
		}
	}

	/**
	 * 得到单元坐标信息
	 * 
	 * @param ele   单元结点编号
	 * @param coord 每个单元结点坐标, 行数为nPerElement, 列数为nDim
	 */
	public void getCoordPerElement(int[] ele, double[][] coord) {
		for (int j = 0; j < nPerElement; j++) {
			coord[j] = nodes[ele[j]].clone();
		}
	}

	/**
	 * 得到边界坐标信息
	 * 
	 * @param ele   边界结点编号
	 * @param coord 每个单元结点坐标, 行数为nPerBoundary, 列数为nDim
	 */
	public void getCoordPerBoundary(int[] ele, double[][] coord) {
		for (int j = 0; j < nPerBoundary; j++) {
			coord[j] = nodes[ele[j]].clone();
		}
	}

	/**
	 * 得到结点坐标
	 * 
	 * @param i 结点编号
	 * @return 结点坐标数组
	 */
	public double[] getNode(int i) {
		return nodes[i].clone();
	}

	public ArrayList<Integer> getNodes(ScalarFunc func, double[] param) {
		ArrayList<Integer> nodeIdx = new ArrayList<>();
		for(int i = 0; i < nv; i++) {
			if(func.action(nodes[i], nodeLabel[i], param) > 0.0) {
				nodeIdx.add(i);
			}
		}
		return nodeIdx;
	}
	/**
	 * 得到单元结点编号
	 * 
	 * @param i 单元编号
	 * @return 单元结点编号数组
	 */
	public int[] getElement(int i) {
		return elements[i].clone();
	}

	public ArrayList<Integer> getElements(ScalarFunc func, double[] param) {
		ArrayList<Integer> eles = new ArrayList<>();
		for (int i = 0; i < nt; i++) {
			int[] ele = elements[i];
			if (func.action(interiorPointInElement(ele), elementLabel[i], param) > 0) {
				eles.add(i);
			}
		}
		return eles;
	}
	public ArrayList<Integer> getElements(int... label) {
		ArrayList<Integer> eles = new ArrayList<>();
		for(int i = 0; i < nt; i++) {
			int eleLabel = getElementLabel(i);
			for(int j = 0, nj = label.length; j < nj; j++) {
				if(eleLabel == label[j]) {
					eles.add(i);
				}
			}
		}
		return eles;
	}

	/**
	 * 得到单元材料编号
	 * 
	 * @param i 单元编号
	 * @return 单元材料编号
	 */
	public int getElementLabel(int i) {
		return elementLabel[i];
	}

	/**
	 * 得到边界结点编号
	 * 
	 * @param i 边界编号
	 * @return 边界编号编号数组
	 */
	public int[] getBoundary(int i) {
		return boundaries[i].clone();
	}

	/**
	 * 得到边界类型编号
	 * 
	 * @param i 边界编号
	 * @return 边界类型编号
	 */
	public int getBoundaryLabel(int i) {
		return boundaryLabel[i];
	}

	/**
	 * 设置网格边界信息
	 * 
	 * @param set HashSet<UnorderIntArray>类, 存储边界结点编号
	 */
	protected void setBoundaries(HashSet<UnorderIntArray> set) {
		this.nb = set.size();
		boundaries = new int[this.nb][];
		boundaryLabel = new int[this.nb];
		int count = 0;
		for (UnorderIntArray ele : set) {
			boundaries[count++] = ele.getIntSet();
		}
	}

	public SIMatrix lineFromElement() {
		SIMatrix edgeMat = new SIMatrix(nv);
		for (int i = 0; i < nt; i++) {
			int[] idx = elements[i];
			for(int j = 0, n = idx.length - 1; j < n; j++) {
				edgeMat.setElement(idx[j], idx[j + 1], i);
			}
			edgeMat.setElement(idx[idx.length - 1], idx[0], i);
		}
		return edgeMat;
	}

	protected void addBoundarys(HashSet<UnorderIntArray> set, int label) {
		this.nb = set.size() + this.nb;
		int[][] newboundarys = new int[this.nb][];
		int[] newboundaryLabel = new int[this.nb];
		int count = 0;
		for (int[] ele : boundaries) {
			newboundarys[count] = ele;
			newboundaryLabel[count] = boundaryLabel[count];
			count++;
		}
		for(UnorderIntArray ele: set) {
			newboundarys[count] = ele.getIntSet();
			newboundaryLabel[count++] = label;
		}
		boundaries = newboundarys;
		boundaryLabel = newboundaryLabel;
	}

	public double[] computeAngles() {
		double[] angles = new double[nv];
		double[][] coordRegion = new double[nPerElement][nDim];
		double[] angle = new double[3];
		for(int i = 0; i < nt; i++) {
			int[] tri = elements[i];
			getCoordPerElement(tri, coordRegion);
			CG.angle2DT3(coordRegion, angle);
			//System.out.println(MVO.toString(coordRegion));
			//System.out.println(MVO.toString(angle));
			for(int j = 0; j < tri.length; j++) {
				angles[tri[j]] += angle[j];
			}
			
		}
		return angles;
	}
	
	/**
	 * 检验网格是否正规
	 * 
	 * @param coord 单元坐标数组
	 * @return 如果计算为正则返回true, 为负则返回false
	 */
	public abstract boolean checkMesh(double[][] coord);

	/**
	 * 对网格作尺度变换
	 * 
	 * @param sxyz 可变实数数组, 要求传入的数都是正数, 长度与网格维数相同
	 */
	public void scale(double... sxyz) {
		int n = sxyz.length;
		ParamCheck.checkEqual(n, nDim);
		for (int i = 0; i < n; i++) {
			ParamCheck.isPositive(sxyz[i]);
		}

		for (int i = 0; i < nv; i++) {
			for (int j = 0; j < nDim; j++) {
				nodes[i][j] *= sxyz[j];
			}
		}
	}

	/**
	 * 对网格作平移变换
	 * 
	 * @param txyz 平衡数组，要求数组长度与网格维数相同
	 */
	public void translate(double... txyz) {
		int n = txyz.length;
		ParamCheck.checkEqual(n, nDim);
		for (int i = 0; i < nv; i++) {
			for (int j = 0; j < nDim; j++) {
				nodes[i][j] += txyz[j];
			}
		}
	}

	/**
	 * 采用传入ScalarFunc函数对网格作变换
	 * 
	 * @param param   变换函数参数
	 * @param funcXYZ 变换函数数组，个数必须与网格维数相同
	 */
	public void transform(double[] param, ScalarFunc... funcXYZ) {
		int n = funcXYZ.length;
		ParamCheck.checkEqual(n, nDim);
		for (int i = 0; i < nv; i++) {
			double[] node = nodes[i].clone();
			for (int j = 0; j < nDim; j++) {
				nodes[i][j] = funcXYZ[j].action(node, 0, param);
			}
		}
	}

	/**
	 * 在给定位移下对网格进行变换，对1维问题，只有第一个数组有效; 对于2维问题，只有前两个数组有效; 对于3维问题，所有三个数组均有效。
	 * 
	 * @param u nv维向量，网格x方向位移
	 * @param v nv维向量，网格y方向位移
	 * @param w nv维向量，网格z方向位移
	 */
	public void deform(double[] u, double[] v, double[] w) {
		switch (nDim) {
		case 1:
			for (int i = 0; i < nv; i++) {
				nodes[i][0] += u[i];
			}
			break;
		case 2:
			for (int i = 0; i < nv; i++) {
				nodes[i][0] += u[i];
				nodes[i][1] += v[i];
			}
			break;
		case 3:
			for (int i = 0; i < nv; i++) {
				nodes[i][0] += u[i];
				nodes[i][1] += v[i];
				nodes[i][2] += w[i];
			}
			break;
		default:
			break;
		}
	}

	/**
	 * 提取指定区域单元编号所对应单元的所有结点编号
	 * 
	 * @param label 单元编号可变数组
	 * @return 对应区域单元编号所在的单元包含的所有结点; 如果没有提供编号，则提取所有结点编号
	 */
	public HashSet<Integer> extractDomainNodes(int... label) {
		HashSet<Integer> set = new HashSet<Integer>();
		if (label.length == 0 || label == null) {
			for (int i = 0; i < nt; i++) {
				int[] ele = elements[i];
				for (int j = 0, m = ele.length; j < m; j++) {
					set.add(ele[j]);
				}
			}
		} else {
			for (int k = 0; k < nt; k++) {
				for (int i = 0, n = label.length; i < n; i++) {
					if (elementLabel[k] == label[i]) {
						int[] ele = elements[k];
						for (int j = 0, m = ele.length; j < m; j++) {
							set.add(ele[j]);
						}
						break;
					}
				}
			}
		}
		return set;
	}

	public Pair<Integer, double[]>[] getNodeCoordinates(int... label) {
		HashSet<Integer> set = extractDomainNodes(label);
		Pair<Integer, double[]>[] result = new Pair[set.size()];
		int idx = 0;
		for(int node: set) {
			result[idx++] = new Pair<>(node, nodes[node]);
		}
		return result;
	}
	/**
	 * 由一个指示函数来提取区域中的某些结点
	 * @param identifyFunc 指示函数，如果区域中结点代入该函数，函数值为正时表示提取该结点，为负则不提取该结点
	 * @param param 指示函数参数
	 * @return HashSet<Integer> 对象, 包含所提取的结点编号集合
	 */
	public HashSet<Integer> extractDomainNodes(ScalarFunc identifyFunc, double[] param) {
		HashSet<Integer> nodesSet = new HashSet<>();
		for(int i = 0; i < nv; i++) {
			double[] coord = nodes[i];
			if(identifyFunc.action(coord, nodeLabel[i], param) > 0) {
				nodesSet.add(i);
			}
		}
		return nodesSet;
	}


//	public SRMatrix renumberingEdges() {
//		SRMatrix A = new SRMatrix(nv);
//		for (int i = 0; i < nt; i++) {
//			int[][] lines = getLinesFromElement(i);
//			for (int j = 0, n = lines.length; j < n; j++) {
//				A.addElement(lines[j][0], lines[j][1], -1.0);
//				A.addElement(lines[j][1], lines[j][0], -1.0);
//			}
//		}
//
//		A.sort();
//		SRMatrix B = A.triU();
//		int nb = B.getNNZ();
//		int[] iIdx = new int[nb];
//		int[] jIdx = new int[nb];
//		double[] values = new double[nb];
//		B.getTriValues(iIdx, jIdx, values);
//		Arrays.setAll(values, i -> i + 1);
//		SRMatrix C = new SRMatrix(iIdx, jIdx, values, nv, nv);
//		return C;
//	}
	/**
	 * 对网格中所有线单元进行重新编号
	 *
	 * @return nv维整型稀疏矩阵, 稀疏矩阵中元素(i, j, value)表示原网格中第i个结点与第j个结点形成的一条边，这条边的编号为value
	 */
	public SIMatrix renumberingEdges() {
		SIMatrix A = new SIMatrix(nv);
		for (int i = 0; i < nt; i++) {
			int[][] lines = getEdgesFromElement(i);
			for (int j = 0, n = lines.length; j < n; j++) {
				A.addElement(lines[j][0], lines[j][1], -1);
				A.addElement(lines[j][1], lines[j][0], -1);
			}
		}

		A.sort();
//		System.out.println(A);
		SIMatrix B = A.triU();
//		System.out.println(B);
		int nb = B.getNNZ();
//		System.out.println(nb);
		int[] iIdx = new int[nb];
		int[] jIdx = new int[nb];
		int[] values = new int[nb];
		B.getTriValues(iIdx, jIdx, values);
		Arrays.setAll(values, i -> i + 1);
		SIMatrix C = new SIMatrix(iIdx, jIdx, values, nv, nv);
		return C;
	}


	public void sortElements() {
		Arrays.sort(elements, (o1, o2)-> {
			//System.out.println("OK");
			double[][] coord1 = new double[nPerElement][nDim];
			double[][] coord2 = new double[nPerElement][nDim];
			getCoordPerElement(o1, coord1);
			getCoordPerElement(o2, coord2);
			double[] p1 = getinteriorPointInElement(coord1);
			double[] p2 = getinteriorPointInElement(coord2);
			//System.out.println(comparators[nDim - 1]);
			return comparators[nDim - 1].compare(p1, p2);
		});
	}

	public Triple<Integer, int[], double[]> pointInMesh(double[] xy, PointInElement pie) {
		this.sortElements();
		int left = 0;
		int right = nt - 1;
		int middle = 0;
		double[][] coord = new double[nPerElement][nDim];
		double[] N = new double[nPerElement];
		int[] ele = null;
		int size = 20;
		while(left <= right) {
			middle = (left + right) / 2;
			ele = elements[middle];
			getCoordPerBoundary(ele, coord);
			double[] p0 = getinteriorPointInElement(coord);
			if (xy[0] < p0[0]) {
				right = middle-1;
			} else {
				left = middle + 1;
			}
			if(right - left < size) {
				break;
			}
		}
//		ele = elements[right];
//		getCoordPerElement(ele, coord);
//		double[] p0 = getinteriorPointInElement(coord);
//		System.out.println(Arrays.toString(ele) + "\t" + Arrays.toString(p0));
//		System.out.println("left:" + left + "\tright: " + right);
		for(int i = left; i < right; i++) {
			ele = elements[i];
			getCoordPerElement(ele, coord);
			double[] p0 = getinteriorPointInElement(coord);
			System.out.println(Arrays.toString(ele) + "\t" + Arrays.toString(p0));
		}
		//System.exit(1);
		for(int i = Math.max(0, left - 20); i < Math.min(right + 20, nt); i++) {
			ele = elements[i];
			getCoordPerElement(ele, coord);
			double[] p0 = getinteriorPointInElement(coord);
			//System.out.println(Arrays.toString(ele) + "\t" + Arrays.toString(p0));
			//System.out.println(MVO.toString(coord));
			//CG.T3ShapeFunction(coord, xy, N);
			//System.out.println(MVO.toString(N));
			if(pie.action(coord, xy, N)) {
				//System.out.println(Arrays.toString(N));
				//System.out.println(Arrays.toString(ele) + "\t" + Arrays.toString(p0));
				return new Triple<Integer, int[], double[]>(i, ele.clone(), N);
			}
		}
		return null;
	}
	
	/**
	 * 由单元编号得到单元每条边条的编号数组，默认为单元编号数组，二维与三维网格需要重写此方法
	 * 
	 * @param i 单元编号
	 * @return 二维数组，得到单元每条边的结点编号
	 */
	public int[][] getEdgesFromElement(int i) {
		return new int[][] { elements[i] };
	}

	/**
	 * 由单元编号得到单元每条边条的编号数组，默认为边界编号数组，一维问题无意义，三维网格需要重写此方法
	 * 
	 * @param i 边界编号
	 * @return 二维数组，得到边界每条边的结点编号
	 */
	public int[][] getLinesFromBoundary(int i) {
		return new int[][] { boundaries[i] };
	}

	/**
	 * 提取指定边界结点集合
	 * 
	 * @param label 网格边界编号可变数组
	 * @return 对应边界编号所在的单元包含的所有结点; 如果没有提供边界编号，则提取所有边界结点编号
	 */
	public HashSet<Integer> extractBoundaryNodes(int... label) {
		HashSet<Integer> set = new HashSet<Integer>();
		if (label == null || label.length == 0) {
			for (int i = 0; i < nb; i++) {
				int[] ele = boundaries[i];
				for (int j = 0, m = ele.length; j < m; j++) {
					set.add(ele[j]);
				}
			}
		} else {
			for (int k = 0; k < nb; k++) {
				for (int i = 0, n = label.length; i < n; i++) {
					if (boundaryLabel[k] == label[i]) {
						int[] idx = boundaries[k];
						for (int j = 0, m = idx.length; j < m; j++) {
							set.add(idx[j]);
						}
					}
				}
			}
		}
		return set;
	}

	/**
	 * 由单元信息提取网格边界集合
	 * 
	 * @return HashSet<UnorderIntArray>对象，存在网格边界信息
	 */
	public HashSet<UnorderIntArray> extractBoundarysFromElements() {
		HashSet<UnorderIntArray> set = new HashSet<>();
		for (int i = 0, n = nt; i < n; i++) {
			int[][] boundarys = getBoundarysFromElement(i);
			for (int j = 0, m = boundarys.length; j < m; j++) {
				UnorderIntArray uia = new UnorderIntArray(boundarys[j]);
				if (!set.contains(uia)) {
					set.add(uia);
				} else {
					set.remove(uia);
				}
			}
		}
		// boundarys.clear();
		// nb = set.size();
		// set.forEach(uia -> boundarys.add(new Line2(uia.getIntSet(), 0)));

		setBoundaries(set);
		return set;
	}


	public HashSet<UnorderIntArray> extractEdgesFromElements() {
		HashSet<UnorderIntArray> set = new HashSet<>();
		for (int i = 0, n = nt; i < n; i++) {
			int[][] boundarys = getBoundarysFromElement(i);
			for (int j = 0, m = boundarys.length; j < m; j++) {
				UnorderIntArray uia = new UnorderIntArray(boundarys[j]);
				set.add(uia);
			}
		}
		return set;
//		UnorderIntArray[] edgesSet = new UnorderIntArray[set.size()];
//		set.toArray(edgesSet);
//		int[][] edges = new int[edgesSet.length][];
//		for(int i = 0; i < edgesSet.length; i++) {
//			edges[i] = edgesSet[i].getIntSet();
//		}
//		return edges;
	}

	/**
	 *
	 * @param domainlabel
	 * @return
	 */
	public HashSet<UnorderIntArray> extractEdgesFromElements(int... domainlabel) {
		HashSet<UnorderIntArray> set = new HashSet<>();
		int nLabel = domainlabel.length;
		for (int i = 0, n = nt; i < n; i++) {
			int label = getElementLabel(i);
			for(int j = 0; j < nLabel; j++) {
				if (domainlabel[j] == label) {
					int[][] boundarys = getBoundarysFromElement(i);
					for (int k = 0, m = boundarys.length; k < m; k++) {
						UnorderIntArray uia = new UnorderIntArray(boundarys[k]);
						set.add(uia);
					}
				}
			}
		}
		return set;
	}

	public HashSet<UnorderIntArray> extractEdgesFromElements(ScalarFunc domainFunc, double[] param) {
		HashSet<UnorderIntArray> set = new HashSet<>();
		for (int i = 0, n = nt; i < n; i++) {
			int[] ele = getElement(i);
			int label = getElementLabel(i);
			double[] xyz = interiorPointInElement(ele);
			if (domainFunc.action(xyz, label, param) > 0.0) {
				int[][] boundarys = getBoundarysFromElement(i);
				for (int k = 0, m = boundarys.length; k < m; k++) {
					UnorderIntArray uia = new UnorderIntArray(boundarys[k]);
					set.add(uia);
				}
			}
		}
		return set;
	}


	public int[][] elementsToEdges() {
		HashSet<UnorderIntArray> set = extractEdgesFromElements();
//		UnorderIntArray[] edgesSet = new UnorderIntArray[set.size()];
//		set.toArray(edgesSet);
//		ArrayList<UnorderIntArray> list = new ArrayList<>();
//		for(int i = 0; i < set.size(); i++) {
//			list.add(edgesSet[i]);
//		}
		ArrayList<UnorderIntArray> list = new ArrayList<>(set);
		int[][] edgeElement = new int[nt][];

		for(int i = 0; i < nt; i++) {
			int[] ele = getElement(i);
			int[][] boundarys = getBoundarysFromElement(i);
			edgeElement[i] = new int[boundarys.length];
			for(int j = 0, nj = boundarys.length; j < nj; j++) {
				edgeElement[i][j] = list.indexOf(new UnorderIntArray(boundarys[j]));
			}
		}
		return edgeElement;
	}

	public void addInterfaceBoundary(ScalarFunc bdFunc, double[] param, int label) {
		HashSet<UnorderIntArray> set = new HashSet<>();
		for (int i = 0, n = nt; i < n; i++) {
			int[][] boundarys = getBoundarysFromElement(i);
			for (int j = 0, m = boundarys.length; j < m; j++) {
				if (bdFunc.action(interiorPointOnBoundary(boundarys[j]) , 0, param) > 0.0) {
					UnorderIntArray uia = new UnorderIntArray(boundarys[j]);
					if (!set.contains(uia)) {
						set.add(uia);
					}
				}
			}
		}
		this.addBoundarys(set, label);
	}

	public int[][] elementsToNode() {
		ArrayList<Integer>[] eleToNode = new ArrayList[nv];
		for(int i = 0; i < nv; i++) {
			eleToNode[i] = new ArrayList<Integer>();
		}
		for(int i = 0; i < nt; i++) {
			int[] ele = getElement(i);
			for(int j = 0, mylen = ele.length; j < mylen; j++) {
				eleToNode[ele[j]].add(i);
			}
		}
		int[][] myEleToNode = new int[nv][];
		for(int i = 0; i < nv; i++) {
			int idx = 0;
			myEleToNode[i] = new int[eleToNode[i].size()];
			for(Integer ele: eleToNode[i]) {
				myEleToNode[i][idx++] = ele;
			}
		}
		return myEleToNode;
	}

	/**
	 * 从部分单元信息提取网格边界集合
	 * 
	 * @return HashSet<UnorderIntArray>对象，存在网格边界信息
	 */
	public void removeElements(ScalarFunc func, double[] param, int newBdLabel) {
		HashSet<UnorderIntArray> bdset = new HashSet<>();
		double[][] coord = new double[nPerElement][nDim];
		boolean[] nodeToRemove = new boolean[nv];
		boolean[] elementToRemove = new boolean[nt];
		int newNtSum = 0;
		for (int i = 0, n = nt; i < n; i++) {
			int[] ele = elements[i];
			getCoordPerElement(ele, coord);
			if(func.action(interiorPointInElement(ele), elementLabel[i], param) < 0.0) {
				int[][] boundarys = getBoundarysFromElement(i);
				for (int j = 0, m = boundarys.length; j < m; j++) {
					UnorderIntArray uia = new UnorderIntArray(boundarys[j]);
					if (!bdset.contains(uia)) {
						bdset.add(uia);
					} else {
						bdset.remove(uia);
					}
				}
				newNtSum++;
			}else {
				elementToRemove[i] = true;
				for(int j = 0; j < nPerElement; j++) {
					nodeToRemove[ele[j]] = true;
				}
			}
		}
		for(int i = 0; i < nt; i++) {
			if(!elementToRemove[i]) {
				for(int j = 0; j < nPerElement; j++) {
					nodeToRemove[elements[i][j]] = false;
				}
			}
		}
		
		int[] old2Newidx = new int[nv];
		int idx = 0;
		for(int i = 0; i < nv; i++) {
			old2Newidx[i] = nodeToRemove[i]?-1:idx++;
		}
		
		int newNv = idx;
		double[][] newNodes = new double[newNv][];
		int[] newNodeLabel = new int[newNv];
		for(int i = 0; i < nv; i++) {
			if(old2Newidx[i] != -1) {
				newNodes[old2Newidx[i]] = nodes[i].clone();
				newNodeLabel[old2Newidx[i]] = nodeLabel[i];
			}
		}
		//System.out.println(MVO.toString(newNodes));
		
		int newNt = newNtSum;
		int[][] newElements = new int[newNt][];
		int[] newElementLabel = new int[newNt];
		int ntIdx = 0;
		int[] newEle = new int[nPerElement];
		for(int i = 0; i < nt; i++) {
			int[] oldEle = elements[i];
			if(!elementToRemove[i]) {
				for(int j = 0; j < nPerElement; j++) {
					newEle[j] = old2Newidx[oldEle[j]];
				}
				newElements[ntIdx] = newEle.clone();
				newElementLabel[ntIdx++] = elementLabel[i];
			}
		}
		//System.out.println(MVO.toString(newElements));
		
		int newBd = bdset.size();
		int bdIdx = 0;
		int[][] newBoundarys = new int[newBd][];
		int[] newBoundaryLabel = new int[newBd];
		for(int i = 0; i < boundaries.length; i++) {
			UnorderIntArray uia = new UnorderIntArray(boundaries[i]);
			if (bdset.contains(uia)) {
				int[] newBdEle = new int[boundaries[i].length];
				for(int j = 0; j < boundaries[i].length; j++) {
					newBdEle[j] = old2Newidx[boundaries[i][j]];
				}
				newBoundarys[bdIdx] = newBdEle;
				newBoundaryLabel[bdIdx++] = boundaryLabel[i];
				bdset.remove(uia);
			}
		}
		
		for(UnorderIntArray array: bdset) {
			int[] oldBdEle = array.getIntSet();
			int[] newBdEle = new int[oldBdEle.length];
			for(int i = 0; i < oldBdEle.length; i++) {
				newBdEle[i] = old2Newidx[oldBdEle[i]];
			}
			newBoundarys[bdIdx] = newBdEle;
			newBoundaryLabel[bdIdx++] = newBdLabel;
		}
		//System.out.println(MVO.toString(newBoundarys));
		
		//System.out.println("New Nv, Nt, Nb: " + newNv + "\t" + newNt + "\t" + newBd);
		
		this.nv = newNv;
		this.nt = newNt;
		this.nb = newBd;
		this.nodes = newNodes;
		this.nodeLabel = newNodeLabel;
		this.elements = newElements;
		this.elementLabel = newElementLabel;
		this.boundaries = newBoundarys;
		this.boundaryLabel = newBoundaryLabel;
	}

	public void removeElements(int[] eles, int newBdLabel) {
		HashSet<UnorderIntArray> bdset = new HashSet<>();
		double[][] coord = new double[nPerElement][nDim];
		boolean[] nodeToRemove = new boolean[nv];
		boolean[] elementToRemove = new boolean[nt];
		int newNtSum = 0;

		for(int i = 0; i < eles.length; i++) {
			int[] ele = getElement(eles[i]);
			elementToRemove[eles[i]] = true;
			for(int j = 0; j < nPerElement; j++) {
				nodeToRemove[ele[j]] = true;
			}
		}


		for(int i = 0; i < nt; i++) {
			if(!elementToRemove[i]) {
				int[][] boundarys = getBoundarysFromElement(i);
				for (int j = 0, m = boundarys.length; j < m; j++) {
					UnorderIntArray uia = new UnorderIntArray(boundarys[j]);
					if (!bdset.contains(uia)) {
						bdset.add(uia);
					} else {
						bdset.remove(uia);
					}
				}
				newNtSum++;
			}
		}

		for(int i = 0; i < nt; i++) {
			if(!elementToRemove[i]) {
				for(int j = 0; j < nPerElement; j++) {
					nodeToRemove[elements[i][j]] = false;
				}
			}
		}

		int[] old2Newidx = new int[nv];
		int idx = 0;
		for(int i = 0; i < nv; i++) {
			old2Newidx[i] = nodeToRemove[i]?-1:idx++;
		}

		int newNv = idx;
		double[][] newNodes = new double[newNv][];
		int[] newNodeLabel = new int[newNv];
		for(int i = 0; i < nv; i++) {
			if(old2Newidx[i] != -1) {
				newNodes[old2Newidx[i]] = nodes[i].clone();
				newNodeLabel[old2Newidx[i]] = nodeLabel[i];
			}
		}
		//System.out.println(MVO.toString(newNodes));

		int newNt = newNtSum;
		int[][] newElements = new int[newNt][];
		int[] newElementLabel = new int[newNt];
		int ntIdx = 0;
		int[] newEle = new int[nPerElement];
		for(int i = 0; i < nt; i++) {
			int[] oldEle = elements[i];
			if(!elementToRemove[i]) {
				for(int j = 0; j < nPerElement; j++) {
					newEle[j] = old2Newidx[oldEle[j]];
				}
				newElements[ntIdx] = newEle.clone();
				newElementLabel[ntIdx++] = elementLabel[i];
			}
		}
		//System.out.println(MVO.toString(newElements));

		int newBd = bdset.size();
		int bdIdx = 0;
		int[][] newBoundarys = new int[newBd][];
		int[] newBoundaryLabel = new int[newBd];
		for(int i = 0; i < boundaries.length; i++) {
			UnorderIntArray uia = new UnorderIntArray(boundaries[i]);
			if (bdset.contains(uia)) {
				int[] newBdEle = new int[boundaries[i].length];
				for(int j = 0; j < boundaries[i].length; j++) {
					newBdEle[j] = old2Newidx[boundaries[i][j]];
				}
				newBoundarys[bdIdx] = newBdEle;
				newBoundaryLabel[bdIdx++] = boundaryLabel[i];
				bdset.remove(uia);
			}
		}

		for(UnorderIntArray array: bdset) {
			int[] oldBdEle = array.getIntSet();
			int[] newBdEle = new int[oldBdEle.length];
			for(int i = 0; i < oldBdEle.length; i++) {
				newBdEle[i] = old2Newidx[oldBdEle[i]];
			}
			newBoundarys[bdIdx] = newBdEle;
			newBoundaryLabel[bdIdx++] = newBdLabel;
		}
		this.nv = newNv;
		this.nt = newNt;
		this.nb = newBd;
		this.nodes = newNodes;
		this.nodeLabel = newNodeLabel;
		this.elements = newElements;
		this.elementLabel = newElementLabel;
		this.boundaries = newBoundarys;
		this.boundaryLabel = newBoundaryLabel;
	}


	public void removeMesh(Mesh myMesh, int newBdLabel) {
		HashSet<UnorderIntArray> bdset = new HashSet<>();
		double[][] coord = new double[nPerElement][nDim];
		double[][] coord2 = new double[nPerElement][nDim];
		boolean[] nodeToRemove = new boolean[nv];
		boolean[] elementToRemove = new boolean[nt];
		ArrayList<Integer> eles = new ArrayList<>();
		int newNtSum = 0;
		int myMeshNt = myMesh.getNt();
		for(int i = 0; i < nt; i++) {
			int[] ele = getElement(i);
			getCoordPerElement(ele, coord);
			double[] p = getinteriorPointInElement(coord);
			for (int j = 0; j < myMeshNt; j++) {
				int[] MyEle = myMesh.getElement(j);
				myMesh.getCoordPerElement(MyEle, coord2);
				double[] myP = myMesh.getinteriorPointInElement(coord2);
				if (MVO.distance(p, myP) < 1.0e-6) {
					eles.add(i);
					break;
				}
			}
		}
		int[] myEles = new int[eles.size()];
		Arrays.setAll(myEles, i -> eles.get(i));
		removeElements(myEles, newBdLabel);
	}

	/**
	 * 将网格导出至文件，可用Tecplot软件打开查看
	 * 
	 * @param fileName 导出文件名
	 */
	public void toTecplot(String fileName) {
		FileWriter fw = null;
		try {
			fw = new FileWriter(fileName);
			fw.write(Tecplot.meshFileString + "\n");
			fw.write(Tecplot.dimVariableString[nDim - 1] + "\n");
			fw.write(Tecplot.ZoneString(this) + "\n");
			for (int i = 0; i < nv; i++) {
				for (int j = 0; j < nDim; j++) {
					fw.write(String.format("%.8f", nodes[i][j]) + "\t ");
				}
				fw.write("\n");
			}

			for (int i = 0; i < nt; i++) {
				for (int j = 0, n = elements[i].length; j < n; j++) {
					fw.write((elements[i][j] + 1) + "\t ");
				}
				fw.write("\n");
			}
			fw.flush();
			fw.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	public void toTecplot(String fileName, double time) {
		if (Math.abs(time) < Constants.Er) {
			try (FileWriter fw = new FileWriter(fileName)) {
				fw.write(Tecplot.solutionFileString + "\n");
				fw.write(Tecplot.dimVariableString[nDim - 1] + "\n");
				fw.write(Tecplot.ZoneStringTime(this) + time + "\n");

				for (int i = 0; i < nv; i++) {
					for (int j = 0; j < nDim; j++) {
						fw.write(nodes[i][j] + "\t");
					}
					fw.write("\n");
				}

				for (int i = 0; i < nt; i++) {
					for (int j = 0, n = elements[i].length; j < n; j++) {
						fw.write((elements[i][j] + 1) + "\t");
					}
					fw.write("\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		} else {
			try (FileWriter fw = new FileWriter(fileName, true)) {
				fw.write(Tecplot.ZoneStringDeformShareTime(this) + time + "\n");
				for (int i = 0; i < nv; i++) {
					for (int j = 0; j < nDim; j++) {
						fw.write((nodes[i][j]) + "\t");
					}
					fw.write("\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		}
	}
	/**
	 * 将网格与计算结果导出至文件，文件可用Tecplot软件打开查看
	 * 
	 * @param fileName 导出文件名
	 * @param x        网格计算结果，实型数组，长度与网格结点数相同
	 */
	public void toTecplot(String fileName, double[] x) {
		FileWriter fw = null;
		try {
			fw = new FileWriter(fileName);
			fw.write(Tecplot.solutionFileString + "\n");
			fw.write(Tecplot.dimVariableString[nDim - 1] + ", u\n");
			fw.write(Tecplot.ZoneString(this) + "\n");
			for (int i = 0; i < nv; i++) {
				for (int j = 0; j < nDim; j++) {
					fw.write(nodes[i][j] + "\t");
				}
				fw.write(x[i] + "\n");
			}

			for (int i = 0; i < nt; i++) {
				for (int j = 0, n = elements[i].length; j < n; j++) {
					fw.write((elements[i][j] + 1) + "\t");
				}
				fw.write("\n");
			}
			fw.flush();
			fw.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}
	
	public void toTecplot(String fileName, Complex[] x) {
		FileWriter fw = null;
		try {
			fw = new FileWriter(fileName);
			fw.write(Tecplot.solutionFileString + "\n");
			fw.write(Tecplot.dimVariableString[nDim - 1] + ", u, v\n");
			fw.write(Tecplot.ZoneString(this) + "\n");
			for (int i = 0; i < nv; i++) {
				for (int j = 0; j < nDim; j++) {
					fw.write(nodes[i][j] + "\t");
				}
				fw.write(String.format("%16.10f\t%16.10f\n", x[i].real, x[i].image));
			}

			for (int i = 0; i < nt; i++) {
				for (int j = 0, n = elements[i].length; j < n; j++) {
					fw.write((elements[i][j] + 1) + "\t");
				}
				fw.write("\n");
			}
			fw.flush();
			fw.close();
		} catch (IOException e) {
			e.printStackTrace();
		}		
	}
	/**
	 * 将网格与计算结果导出至文件，文件可用Tecplot软件打开查看
	 * 
	 * @param fileName 导出文件名
	 * @param mat        网格计算结果，实型二维数组，行数等于计算结果个数，列数与网格结点数相同
	 */
	public void toTecplot(String fileName, double[][] mat) {
		int col = mat.length;
		try (FileWriter fw = new FileWriter(fileName)) {
			fw.write(Tecplot.solutionFileString + "\n");
			fw.write(Tecplot.dimVariableString[nDim - 1]);
			for (int i = 0; i < col; i++) {
				fw.write(", \"u" + (i + 1) + "\"");
			}
			fw.write("\n");

			fw.write(Tecplot.ZoneString(this) + "\n");
			for (int i = 0; i < nv; i++) {
				for (int j = 0; j < nDim; j++) {
					fw.write(nodes[i][j] + "\t ");
				}
				int n = mat.length;
				for (int j = 0; j < n; j++) {
					fw.write(mat[j][i] + "\t ");
				}
				fw.write("\n");
			}
			for (int i = 0; i < nt; i++) {
				for (int j = 0, n = elements[i].length; j < n; j++) {
					fw.write((elements[i][j] + 1) + "\t ");
				}
				fw.write("\n");
			}
			fw.flush();
			fw.close();

		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	/**
	 * @brief 将各个时刻计算出的解导出到至文件，以供tecplot文件读取
	 * 
	 * 用法: 在进行时间积分前，需要先调用一次该函数保存在time = 0 时刻的解如
	 *     ...
	 *     mesh.toTecplot(fileName, 0.0, x);
	 *     ...
	 *     在时间上积分中每计算出一步在时间time = t上的有限元解x时，调用
	 *     ...
	 *     mesh.toTecplot(fileName, t, x);
	 *     ...
	 *     保存当前时刻的解，注意两处调用的文件名必须相同
	 * 
	 * @param fileName 文件名
	 * @param time 时间变量
	 * @param x 特定时刻时的解向量
	 */
	public void toTecplot(String fileName, double time, double[] x) {
		if (Math.abs(time) < Constants.Er) {
			try (FileWriter fw = new FileWriter(fileName)) {
				fw.write(Tecplot.solutionFileString + "\n");
				fw.write(Tecplot.dimVariableString[nDim - 1] + ", u\n");
				fw.write(Tecplot.ZoneStringTime(this) + time + "\n");
				for (int i = 0; i < nv; i++) {
					for (int j = 0; j < nDim; j++) {
						fw.write(nodes[i][j] + "\t");
					}
					fw.write(x[i] + "\n");
				}

				for (int i = 0; i < nt; i++) {
					for (int j = 0, n = elements[i].length; j < n; j++) {
						fw.write((elements[i][j] + 1) + "\t");
					}
					fw.write("\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		} else {
			try (FileWriter fw = new FileWriter(fileName, true)) {
				fw.write(Tecplot.ZoneStringShareTime(this) + time + "\n");
				for (int i = 0, n = x.length; i < n; i++) {
					fw.write(x[i] + "\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		}
	}

	/**
	 * @brief 将各个时刻计算出的解导出到至文件，以供tecplot文件读取
	 * @param fileName 文件名
	 * @param time 时间变量
	 * @param mat 二维数组, 行数为解的个数, 列数为网格结点数 
	 */
	public void toTecplot(String fileName, double time, double[][] mat) {
		int row = mat.length;
		if (Math.abs(time) < Constants.Er) {
			try (FileWriter fw = new FileWriter(fileName)) {
				fw.write("title = mesh file\n");
				fw.write(Tecplot.dimVariableString[nDim - 1]);
				for (int j = 0; j < row; j++) {
					fw.write(", \"u" + (j + 1) + "\"");
				}
				fw.write("\n");

				fw.write(Tecplot.ZoneStringTime(this) + time + "\n");
				for (int i = 0; i < nv; i++) {
					for (int j = 0; j < nDim; j++) {
						fw.write(nodes[i][j] + "\t ");
					}
					for (int j = 0; j < row; j++) {
						fw.write(mat[j][i] + "\t ");
					}
					fw.write("\n");
				}
				for (int i = 0; i < nt; i++) {
					for (int j = 0, n = elements[i].length; j < n; j++) {
						fw.write((elements[i][j] + 1) + "\t ");
					}
					fw.write("\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		} else {
			try (FileWriter fw = new FileWriter(fileName, true)) {
				fw.write(Tecplot.ZoneStringShareTime(this) + time + "\n");
				int col = mat[0].length;
				for (int j = 0; j < col; j++) {
					for (int i = 0; i < row; i++) {
						fw.write(mat[i][j] + "\t");
					}
					fw.write("\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		}
	}

	
	public void toTecplotDeform(String fileName, double[][] disp, double coef) {
		try (FileWriter fw = new FileWriter(fileName)) {
			fw.write(Tecplot.solutionFileString + "\n");
			fw.write(Tecplot.dimVariableString[nDim - 1] + "\n");
			fw.write(Tecplot.ZoneString(this) + "\n");
			for (int i = 0; i < nv; i++) {
				for (int j = 0; j < nDim; j++) {
					fw.write((nodes[i][j] + disp[j][i] * coef) + "\t");
				}
				fw.write("\n");
			}

			for (int i = 0; i < nt; i++) {
				for (int j = 0, n = elements[i].length; j < n; j++) {
					fw.write((elements[i][j] + 1) + "\t");
				}
				fw.write("\n");
			}
			fw.flush();
			fw.close();
		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	public void toTecplotDeform(String fileName, double time, double[][] disp, double coef) {
		if (Math.abs(time) < Constants.Er) {
			try (FileWriter fw = new FileWriter(fileName)) {
				fw.write(Tecplot.solutionFileString + "\n");
				fw.write(Tecplot.dimVariableString[nDim - 1] + "\n");
				fw.write(Tecplot.ZoneStringTime(this) + time + "\n");

				for (int i = 0; i < nv; i++) {
					for (int j = 0; j < nDim; j++) {
						fw.write(nodes[i][j] + "\t");
					}
					fw.write("\n");
				}

				for (int i = 0; i < nt; i++) {
					for (int j = 0, n = elements[i].length; j < n; j++) {
						fw.write((elements[i][j] + 1) + "\t");
					}
					fw.write("\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		} else {
			ParamCheck.checkEqual(disp.length, nDim);
			try (FileWriter fw = new FileWriter(fileName, true)) {
				fw.write(Tecplot.ZoneStringDeformShareTime(this) + time + "\n");
				for (int i = 0; i < nv; i++) {
					for (int j = 0; j < nDim; j++) {
						fw.write((nodes[i][j] + disp[j][i] * coef) + "\t");
					}
					fw.write("\n");
				}
				fw.flush();
				fw.close();
			} catch (Exception e) {
				e.printStackTrace();
			}
		}
	}

	void toGmshSparse(String fileName, double[] u) {
		assert(u.length == nv);
		try (FileWriter fw = new FileWriter(fileName)) {
			fw.write("View \"A scalar field\" {\n");
			double[][] coord = new double[nPerElement][nDim];
			for(int i = 0; i < nt; i++) {
				fw.write(gmshType[0] + "(");
				int[] ele = getElement(i);
				getCoordPerElement(ele, coord);
				for(int j = 0; j < nPerElement; j++) {
					if (nDim == 3) {
						fw.write(coord[j][0] + ", "+coord[j][0] + ", "+coord[j][0] + "");
					}
				}
			}
			fw.write("}\n");
		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	public void toGmsh(String fileName) {
		try(FileWriter fw = new FileWriter(fileName)) {
			fw.write("$MeshFormat\n");
			fw.write("2.2 0 8\n");
			fw.write("$EndMeshFormat\n");
			fw.write("$Nodes\n");
			fw.write( nv +"\n");
			for(int i = 0; i < nv; i++) {
				fw.write((i + 1) + "  ");
				for(int j = 0; j < nDim; j++) {
					fw.write(nodes[i][j] + "  ");
				}
				for(int j = nDim; j < 3; j++) {
					fw.write("0.0 ");
				}
				fw.write("\n");
			}
			fw.write("$EndNodes\n");
			fw.write("$Elements\n");
			fw.write(nt + "\n");
			for(int i = 0; i < nt; i++) {
				fw.write((i + 1) + "  " + gmshElementType + "  2 " + elementLabel[i] + " 2 ");
				for(int j = 0; j < nPerElement; j++) {
					fw.write((elements[i][j] + 1) + "  ");
				}
				fw.write("\n");
			}
			fw.write("$EndElements\n");
		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	public void toGmsh(String fileName, double[] u) {
		try(FileWriter fw = new FileWriter(fileName)) {
			fw.write("$MeshFormat\n");
			fw.write("2.2 0 8\n");
			fw.write("$EndMeshFormat\n");
			fw.write("$Nodes\n");
			fw.write( nv +"\n");
			for(int i = 0; i < nv; i++) {
				fw.write((i + 1) + "  ");
				for(int j = 0; j < nDim; j++) {
					fw.write(nodes[i][j] + "  ");
				}
				for(int j = nDim; j < 3; j++) {
					fw.write("0.0 ");
				}
				fw.write("\n");
			}
			fw.write("$EndNodes\n");
			fw.write("$Elements\n");
			fw.write(nt + "\n");
			for(int i = 0; i < nt; i++) {
				fw.write((i + 1) + "  " + gmshElementType + "  2 " + elementLabel[i] + " 2 ");
				for(int j = 0; j < nPerElement; j++) {
					fw.write((elements[i][j] + 1) + "  ");
				}
				fw.write("\n");
			}
			fw.write("$EndElements\n");
			fw.write("$NodeData\n");
			fw.write("1\n");
			fw.write("Scalar view\n");
			fw.write("1\n0.0\n3\n0\n1\n");
			fw.write(nv + "\n");
			for(int i = 0; i < nv; i++) {
				fw.write((i + 1) + "  " + u[i] + "\n");
			}
			fw.write("$EndNodeData\n");
		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	public void CMKRelabelling() {
		SIMatrix graph = constructAGraph();
		boolean[] flag = new boolean[nv];
		int[] degree = new int[nv];
		Arrays.setAll(degree, i -> graph.mat[i].size());
		Pair<Integer, Integer> pair = new Pair<>();
		MVO.min(degree, pair);
		Vector<Integer> queue = new Vector<>();
		int startNode = pair.getSecond();
		flag[startNode] = true;
		queue.add(startNode);
		Vector<Pair<Integer, Integer>> w = new Vector<>();
		int head = 0;
		while (queue.size() < nv) {
			for (SIEntry ele : graph.mat[startNode]) {
				if (!flag[ele.col]) {
					w.add(new Pair<Integer, Integer>(ele.col, degree[ele.col]));
					flag[ele.col] = true;
				}
			}
			w.sort((a, b) -> a.second > b.second ? 1 : -1);
			for (Pair<Integer, Integer> e : w) {
				queue.add(e.getFirst());
			}
			w.clear();
			startNode = queue.get(++head);
		}

		Collections.reverse(queue);

		double[][] newNodes = new double[nv][nDim];
		for (int i = 0; i < nv; i++) {
			newNodes[i] = nodes[queue.get(i)].clone();
		}

		nodes = newNodes;

		for (int i = 0; i < nt; i++) {
			int[] idx = elements[i];
			int lenIdx = idx.length;
			for (int j = 0; j < lenIdx; j++) {
				idx[j] = queue.indexOf(idx[j]);
			}
		}

		for (int i = 0; i < nb; i++) {
			int[] idx = boundaries[i];
			int lenIdx = idx.length;
			for (int j = 0; j < lenIdx; j++) {
				idx[j] = queue.indexOf(idx[j]);
			}
		}
		// queue.stream().forEach(System.out::println);
	}

	public SIMatrix constructAGraph() {
		SIMatrix result = new SIMatrix(nv);
		for (int k = 0; k < nt; k++) {
			int[] ele = elements[k];
			int n = ele.length;
			for (int i = 0; i < n; i++) {
				for (int j = 0; j < n; j++) {
					if (i != j) {
						result.setElement(ele[i], ele[j], 1);
					}
				}
			}
		}
		return result;
	}

	/**
	 * 
	 * 将原网格中加上已有的网格，产生新的网格, 一般要求两网格类型相同, 网格边界能相互匹配
	 * 
	 * @param meshB
	 *            网格B
	 */
	public void add(Mesh meshB) {
		double tol = 1.0e-3;
		boolean flag = false;
		int nvA, nvB, nConnected, n, nBd;
		int[] newNodesB = new int[meshB.nv];
		int[] newBdB = new int[meshB.nb];
		int[] ABd = new int[nb];
		Arrays.fill(newNodesB, -1);
		Arrays.fill(newBdB, -1);
		Arrays.fill(ABd, -1);
		nvA = nv;
		nvB = meshB.nv;
		for (int i = 0; i < nvA; i++) {
			double[] xyzA = nodes[i];
			for (int j = 0; j < nvB; j++) {
				double[] xyzB = meshB.nodes[j];
				if (MVO.distance(xyzA, xyzB) < tol) {
					newNodesB[j] = i;
					break;
				}
			}
		}

		for (int i = 0; i < meshB.nb; i++) {
			int[] idx = meshB.boundaries[i];
			int[] temp = new int[idx.length];
			Arrays.setAll(temp, ii -> newNodesB[idx[ii]]);
			UnorderIntArray uia1 = new UnorderIntArray(temp);
			flag = true;
			for (int j = 0, tn = temp.length; j < tn; j++) {
				if (temp[j] <= -1) {
					flag = false;
				}
			}
			if (flag) {
				newBdB[i] = 1;
				for (int j = 0; j < nb; j++) {
					int[] idx2 = boundaries[j];
					UnorderIntArray uia2 = new UnorderIntArray(idx2);
					if (uia1.equals(uia2)) {
						ABd[j] = 1;
						break;
					}
				}
			}
		}

		// nConnected = (int)(IntStream.of(newNodesB).filter(iidx -> iidx != -1).count());
		// nBd = (int)(IntStream.of(newBdB).filter(iidx -> iidx != -1).count());
		nConnected = 0;
		n = nvA;
		for (int i = 0; i < nvB; i++) {
			// 如果B中的结点没有与网格中的结点配对，则设置网格B中的结点在新的网格中的结点编号
			if (newNodesB[i] == -1) {
				newNodesB[i] = n;
				n++;
			} else {
				// 如果已经配对，则计算相配对的结点的个数
				nConnected++;
			}

		}

		nBd = 0;
		n = nb;
		for (int i = 0; i < meshB.nb; i++) {
			if (newBdB[i] == -1) {
				newBdB[i] = n;
				n++;
			} else {
				nBd++;
			}
		}

		int i = 0;
		int ntA = nt;
		int nbA = nb;
		nv = nvA + nvB - nConnected;
		nt = nt + meshB.nt;
		nb = nb + meshB.nb - nBd;
		double[][] newNodes = new double[nv][];
		int[] newNodeLabel = new int[nv];
		int[][] newElements = new int[nt][];
		int[] newElementLabel = new int[nt];
		int[][] newBoundarys = new int[nb][];
		int[] newBoundaryLabel = new int[nb];
		for (i = 0; i < nvA; i++) {
			newNodes[i] = this.nodes[i];
			newNodeLabel[i] = this.nodeLabel[i];
		}

		int count = nvA;
		for (i = 0; i < nvB; i++) {
			if (newNodesB[i] > nvA - 1) {
				double[] node = meshB.nodes[i];
				newNodes[count] = node;
				newNodeLabel[count++] = meshB.nodeLabel[i];
			}
		}

		for (i = 0; i < ntA; i++) {
			newElements[i] = this.elements[i];
			newElementLabel[i] = this.elementLabel[i];
		}

		count = ntA;
		for (i = 0; i < meshB.nt; i++) {
			int[] ele = meshB.elements[i];
			int[] newEle = new int[ele.length];
			Arrays.setAll(newEle, ii -> newNodesB[ele[ii]]);
			newElements[count] = newEle;
			newElementLabel[count++] = meshB.elementLabel[i];
		}

		for (i = 0; i < nbA; i++) {
			if (ABd[i] == 1) {
				boundaryLabel[i] = 0;
			}
			newBoundarys[i] = boundaries[i];
			newBoundaryLabel[i] = boundaryLabel[i];
		}
		// System.out.println(Arrays.toString(newBdB));

		count = nbA;
		for (i = 0; i < meshB.nb; i++) {
			if (newBdB[i] != 1) {
				int[] li = meshB.boundaries[i];
				int[] newEle = new int[li.length];
				Arrays.setAll(newEle, ii -> newNodesB[li[ii]]);

				newBoundarys[count] = newEle;
				newBoundaryLabel[count++] = meshB.boundaryLabel[i];
			}
		}

		this.nodes = newNodes;
		this.nodeLabel = newNodeLabel;
		this.elements = newElements;
		this.elementLabel = newElementLabel;
		this.boundaries = newBoundarys;
		this.boundaryLabel = newBoundaryLabel;
	}

	/**
	 * 将多个网格一起加到原网格中
	 * @param mesh 网格数组
	 */
	public void add(Mesh... mesh) {
		for (int i = 0; i < mesh.length; i++) {
			this.add(mesh[i]);
		}
	}
	
	public int[][] extractPeriodicBoundary(int label1, int label2, Direct direct) {
		HashSet<Integer> bd1 = this.extractBoundaryNodes(label1);
		//System.out.println(bd1.toString());
		HashSet<Integer> bd2 = this.extractBoundaryNodes(label2);
		//System.out.println(bd2.toString());
		int n = bd1.size();
		ParamCheck.checkEqual(n, bd2.size());
		int[][] idxPair = new int[n][2];
		int count = 0;
		double[] coord1 = new double[nDim];
		double[] coord2 = new double[nDim];
		int directValue = direct.getValue();
		for(int bdNode1:bd1) {
			idxPair[count][0] = bdNode1;
			this.getCoordPerNode(bdNode1, coord1);
			double[] reducedCoord1 = MVO.removeIdx(coord1, directValue);
			for(int bdNode2:bd2) {
				this.getCoordPerNode(bdNode2, coord2);
				double[] reducedCoord2 = MVO.removeIdx(coord2, directValue);
				if(MVO.distance(reducedCoord1, reducedCoord2) < 1.0e-6) {
					idxPair[count++][1] = bdNode2;
					bd2.remove(bdNode2);
					break;
				}
			}
		}
		//System.out.println("Count = " + count);
		return idxPair;
	}

	public int[][] extractPeriodicBoundary(int label1, int label2, double[] tangent) {
		HashSet<Integer> bd1 = this.extractBoundaryNodes(label1);
		HashSet<Integer> bd2 = this.extractBoundaryNodes(label2);
		int n = bd1.size();
		ParamCheck.checkEqual(n, bd2.size());
		int[][] idxPair = new int[n][2];
		int count = 0;
		double[] coord1 = new double[nDim];
		double[] coord2 = new double[nDim];
		for(int bdNode1:bd1) {
			idxPair[count][0] = bdNode1;
			this.getCoordPerNode(bdNode1, coord1);
			//System.out.println("coord1:");
			//System.out.println(Arrays.toString(coord1));
			//System.out.println("coord2:");
			for(int bdNode2:bd2) {
				this.getCoordPerNode(bdNode2, coord2);
				//System.out.println(Arrays.toString(coord2));
				double[] normal = MVO.add(coord1, -1.0, coord2);
				//System.out.println("normal: " + Arrays.toString(normal));
				//System.out.println(MVO.dot_product(tangent, normal));
				if(Math.abs(MVO.dot_product(tangent, normal)) < 1.0e-8) {
					idxPair[count++][1] = bdNode2;
					bd2.remove(bdNode2);
					break;
				}
			}
			//System.exit(1);
		}
		return idxPair;
	}
	
	public int[][] extractPeriodicBoundary(int[][] label, Direct[] direct) {
		int len = label.length;
		int cc = 0;
		ArrayList<int[][]> pairs = new ArrayList<>();
		for(int i = 0; i < len; i++) {
			HashSet<Integer> bd1 = this.extractBoundaryNodes(label[i][0]);
			HashSet<Integer> bd2 = this.extractBoundaryNodes(label[i][1]);
			int n = bd1.size();
			ParamCheck.checkEqual(n, bd2.size());
			int[][] idxPair = new int[n][2];
			cc += n;
			int count = 0;
			double[] coord1 = new double[nDim];
			double[] coord2 = new double[nDim];
			int directValue = direct[i].getValue();
			for(int bdNode1:bd1) {
				idxPair[count][0] = bdNode1;
				this.getCoordPerNode(bdNode1, coord1);
				double[] reducedCoord1 = MVO.removeIdx(coord1, directValue);
				for(int bdNode2:bd2) {
					this.getCoordPerNode(bdNode2, coord2);
					double[] reducedCoord2 = MVO.removeIdx(coord2, directValue);
					if(MVO.distance(reducedCoord1, reducedCoord2) < Constants.Er) {
						idxPair[count++][1] = bdNode2;
						bd2.remove(bdNode2);
						break;
					}
				}
			}
			pairs.add(idxPair);
		}
		int[][] result = new int[cc][];
		int idx = 0;
		for(int[][] ele: pairs) {
			for(int i = 0, n = ele.length; i < n; i++) {
				result[idx++] = ele[i];
			}
		}
		return result;
	}

	/**
	 * 使用局部协调RML(RC-PML)方法将PML区域结点坐标设置为其复数形式
	 * @param PMLDomainLabel PML区域编号
	 * @param k 波数
	 * @param inLabel PML区域内边界编号
	 * @param outLabel PML区域外边界编号
	 * @return HashMap对象, 键为整数, 表示结点编写, 值为复数, 表示结点对应的复数坐标形式
	 */
	public HashMap<Integer, Complex[]> setPMLCoordinates(double k, int[] inLabel, int[] outLabel, int... PMLDomainLabel) {
		Pair<Integer, double[]>[] PMLnodes = getNodeCoordinates(PMLDomainLabel);
		HashSet<Integer> inBdNodes = extractBoundaryNodes(inLabel);
		double[][] inXY = new double[inBdNodes.size()][];
		int idx = 0;
		for(int node: inBdNodes) {
			inXY[idx++] = nodes[node];
		}
		HashSet<Integer> outBdNodes = extractBoundaryNodes(outLabel);
		double[][] outXY = new double[outBdNodes.size()][];
		idx = 0;
		for(int node:outBdNodes) {
			outXY[idx++] = nodes[node];
		}
		HashMap<Integer, Complex[]> myMap = new HashMap<>();
		for(int i = 0; i < PMLnodes.length; i++) {
			Complex[] newXY = LC_PML(PMLnodes[i].getSecond(), k, inXY, outXY);
			myMap.put(PMLnodes[i].getFirst(), newXY);
		}
		return myMap;
	}

	/**
	 * 计算PML区域结点的复数形式
	 * @param xy PML区域结点坐标值
	 * @param k 波数
	 * @param inXY PML区域内边界结点坐标数组
	 * @param outXY PML区域外边界结点坐标数组
	 * @return xy坐标对应的复数形式
	 */
	public Complex[] LC_PML(double[] xy, double k, double[][] inXY, double[][] outXY) {
		Complex[] complexXY = new Complex[this.nDim];
		//参数alpha值
		double alpha = 7.0 * k;
		//PML衰减率, 整数2或3
		int m = 3;
		//找到在内PML边界最接近点P的点, 定义为xy0
		int inLength = inXY.length;
		double[] dis = new double[inLength];
		Arrays.setAll(dis, i-> MVO.distance(inXY[i], xy));
		Pair<Integer, Double> myPair = new Pair<>();
		MVO.min(dis, myPair);
		double[] xy0 = inXY[myPair.getFirst()];

		//找到外PML边界上沿单位向量方向上的点
		double ksi = myPair.getSecond();
		int outLength = outXY.length;
		//vpXY: 从xy0到xy1的单位向量
		double[][] vpXY = new double[outLength][];
		for(int i = 0; i < outLength; i++) {
			vpXY[i] = MVO.add(outXY[i], -1.0, xy0);
			double norm = MVO.L2Norm(vpXY[i]);
			for(int j = 0; j < this.nDim; j++) {
				vpXY[i][j] /= norm;
			}
		}

		double[] vxy = new double[this.nDim];
		Arrays.setAll(vxy, i -> xy[i] - xy0[i]);
		double norm_vxy = MVO.L2Norm(vxy);
		//nxy: 从xy0到xy的单位向量
		double[] nxy = new double[this.nDim];
		Arrays.setAll(nxy, i-> vxy[i] / norm_vxy);
		if (norm_vxy < 1.e-8) {
			Arrays.setAll(complexXY, i->new Complex(xy[i], 0.0));
		} else {
			//找到最接近的那个点xy1
			double[] temp = new double[outLength];
			Arrays.setAll(temp, i -> Math.acos(MVO.dot_product(vpXY[i], nxy)));
			Pair<Integer, Double> myPair2 = new Pair<>();
			MVO.min(temp, myPair2);
			double[] xy1 = outXY[myPair2.getFirst()];
			//局部pml厚度
			double dpml = MVO.distance(xy0, xy1);
			for(int i = 0; i < this.nDim; i++) {
				complexXY[i] = new Complex(xy[i], -alpha/k * Math.pow(ksi, m) /(m * Math.pow(dpml, m - 1)) * nxy[i]);
			}
		}
		return complexXY;
	}

	public HashSet<Integer>[] SuperconvergencePatches(int minimumPatch) {
		HashSet<Integer>[] node2Elements = new HashSet[nv];
		for(int i = 0; i < nv; i++) {
			node2Elements[i] = new HashSet<Integer>();
		}
		for(int i = 0; i < nt; i++) {
			int[] ele = getElement(i);
			for(int idx : ele) {
				node2Elements[idx].add(i);
			}
		}

//		for(int i = 0; i < nv; i++) {
//			System.out.println(node2Elements[i]);
//		}
		HashSet<Integer>[] newNode2Elements = new HashSet[nv];
		for(int i = 0; i < nv; i++) {
			newNode2Elements[i] = new HashSet<Integer>();
		}
		for(int i = 0; i < nv; i++) {
			for(int idx : node2Elements[i]) {
				newNode2Elements[i].add(idx);
			}
			while(newNode2Elements[i].size() < minimumPatch) {
				Integer[] thisIdx = new Integer[newNode2Elements[i].size()];
				node2Elements[i].toArray(thisIdx);
				for(int iidx: thisIdx) {
					int[] ele = getElement(iidx);
					for(int j = 0, n = ele.length; j < n; j++) {
						for(int otherIdx : node2Elements[ele[j]]) {
							newNode2Elements[i].add(otherIdx);
						}
					}
				}
			}
		}
		return newNode2Elements;
	}
	/**
	 * 网格结点数
	 */
	protected int nv = 0;

	/**
	 * 网格单元数
	 */
	protected int nt = 0;

	/**
	 * 网格边界数
	 */
	protected int nb = 0;

	/**
	 * 网格维数
	 */
	protected int nDim = 0;

	/**
	 * 每个单元结点个数
	 */
	protected int nPerElement = 0;

	/**
	 * 每个边界结点个数
	 */
	protected int nPerBoundary = 0;

	/**
	 * 每个单元中的边界数
	 */
	protected int nBoundaryPerElement = 0;

	/**
	 * 每个单元中的一维边数
	 */
	protected int nEdgePerElement = 0;

	/**
	 * 每个边界中的一维边数
	 */
	protected int nEdgePerBoundary = 0;

	/**
	 * type of element in the software "tecplot", used for post-processing
	 */
	protected String tecplotType = null;
	protected String[] meditMeshType = null;
	protected String[] gmshType = null;
	protected int gmshElementType;

	/**
	 * 网格结点数组，维数: (nv, ndim)
	 */
	double[][] nodes;

	/**
	 * 网格结点编号, 长度: nv
	 */
	int[] nodeLabel;

	/**
	 * 网格单元编号数组, 维数: (nt, nPerElement)
	 */
	int[][] elements;

	/**
	 * 网格单元材料编号, 长度: nt
	 */
	int[] elementLabel;

	/**
	 * 网格边界编号数组, 维数: (nb, nPerBoundary)
	 */
	int[][] boundaries;

	/**
	 * 网格边界编号, 长度: nb
	 */
	int[] boundaryLabel;
	
	private Comparator<double[]> comparator1 = (o1, o2) -> o1[0]-o2[0]>0?1:-1;
	private Comparator<double[]> comparator2 = (o1, o2) -> {
			if (Math.abs(o1[0] - o2[0]) < Constants.Er) {
				return o1[1] - o2[1] > 0?1:-1;
			}
			return o1[0] - o2[0] > 0?1:-1;};
	private Comparator<double[]> comparator3 = (p1, p2) -> {
		if (Math.abs(p1[0] - p2[0]) < Constants.Er) {
			if (Math.abs(p1[1] - p2[1]) < Constants.Er) {
				return p1[2] - p2[2] > 0?1:-1;
			}
			return p1[1] - p2[1] > 0?1:-1;
		}
		return p1[0] - p2[0] > 0?1:-1;
	};
	private Comparator<double[]>[] comparators;



}
